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We investigate numerically simulated collisions between experimentally realistic Bose-Einstein condensate 
wavepackets, within a regime where highly populated scattering haloes are formed. The theoretical basis for 
this work is the truncated Wigner method, for which we present a detailed derivation, paying particular attention 
to its validity regime for colliding condensates. This paper is an extension of our previous Letter [7], and we 
investigate both single-trajectory solutions, which reveal the presence of quantum turbulence in the scattering 
halo, and ensembles of trajectories, which we use to calculate quantum-mechanical correlation functions of the 
field. 
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I. INTRODUCTION 

In the same way as a classical electromagnetic field obey- 
ing Maxwell's equations arises as an assembly of photons all 
in the same quantum state, a Bose-Einstein condensate, com- 
posed of Bosonic atoms all in the same quantum state, be- 
haves very much like a classical field ^'(x, t), whose equation 
of motion is the Gross-Pitaevskii equation 
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The Gross-Pitaevskii equation has been extraordinarily suc- 
cessful in describing a wide range of the phenomena associ- 
ated with Bose-Einstein condensates — nevertheless, there are 
phenomena in which the quantized nature of this field is im- 
portant. For example, when two Bose-Einstein condensates 
collide at a sufficiently high velocity, a halo of elastically scat- 
tered atoms is produced [1-3]. The Gross-Pitaevskii equation 
with initial conditions corresponding to two Bose-Einstein 
condensates does not predict this scattering — it is a direct ef- 
fect of the fact that the quantized field consists of interacting 
particles, and at its most elementary level, this halo is simply 
the result of elastic scattering of the constituent particles in the 
Bose-Einstein condensate. 

A description of this phenomenon by means of the phe- 
nomenological inclusion of loss terms [4], in a way remi- 
niscent of the Boltzmann equation, was successful for colli- 
sions of smaller condensates, but seriously under-estimated 
the scattering for higher condensate densities. This underes- 
timation arises because at high condensate densities, the final 
states into which the scattering occurs can become sufficiently 
highly occupied to cause Bosonic stimulation, which a simple 
Boltzmann treatment cannot produce. On the other hand, a 
treatment in terms of linearized quantum field theory [5, 6] 
can deal with the Bosonic stimulation effects of highly occu- 
pied final states, but, being linearized, can treat neither the 
effects of large depletion, nor the essentially classical nonlin- 
earity effects so well handled by the Gross-Pitaevskii equa- 
tion. 

An alternative method is the truncated Wigner method. 
This is a treatment of quantum field theory in which quantum 
mechanics is simulated by a classical random process. The 
method is approximate, but nevertheless very useful, both for 



quantum optical systems, where the field under consideration 
is the electromagnetic field, and for Bose-Einstein conden- 
sates, where the matter-wave field is under consideration, and 
reproduces many quantum mechanical features, such as quan- 
tum correlation functions, of these systems. The application 
to Bose-Einstein condensates has been developed by several 
groups [10, 13-23] with considerable success. 

In qualitative terms, the truncated Wigner method pro- 
vides a description of quantum field theory based on the 
Gross-Pitaevskii equationin which quantum mechanical vac- 
uum fluctuations are simulated by by adding appropriate clas- 
sical fluctuations in addition to the coherent field of the initial 
state of the Bose-Einstein condensate. These amount to half 
a quantum per degree of freedom, corresponding to the zero 
point energy of the harmonic oscillator which represents each 
mode of the field — the precise way in which this is done is 
presented in Sect.IIC2. The elastic scattering effects which 
are not produced directly by a solution of the Gross-Pitaevskii 
equation for a coherent initial condensate field now appear as 
four-wave mixing between the coherent condensate field and 
the simulated vacuum fluctuations. 

Since the number of modes in a local quantum field theory 
is infinite, the addition of a half a quantum per mode does in 
principle introduce a infinite density of vacuum fluctuations. 
Nevertheless, this does not invalidate the method. In practice 
what is used for the description of a Bose-Einstein conden- 
sate is an effective field theory[&, 9], valid only on a rather 
coarse spatial scale. Because of this, one can use a quantum 
field theory which contains only modes up to a certain cutoff 
value of the momentum, with an appropriately renormalized 
interaction. In a practical implementation of a cutoff field the- 
ory, one uses the standard contact interaction, with strength 
proportional to the scattering length, with corrections which 
depend on the size of the momentum cutoff. The density of 
added vacuum fluctuations is then quite finite, but cutoff de- 
pendent. The cutoff dependent effects of the vacuum fluc- 
tuations are then compensated by the cutoff dependent inter- 
action strength. In simulations we must introduce this cut- 
off explicitly — in practice, for the choices of parameters we 
make, the cutoff corrections are 0.5% — see Sects. II A 2, IV A. 

A related approach, the positive-P method [10, 25], (and its 
generalization, the gauge positive-P method) has also proven 
useful in describing Bose-condensed systems. Both the trun- 
cated Wigner and positive-P methods are examples of the gen- 
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eral class of "classical field methods" well known in quantum 
optics [26]. The positive-P methods are in principle exact, but 
are still difficult to implement for the kind of system we are 
considering here. However, technical progress is currently be- 
ing made [11, 12], and there does seem to be real promise that 
these exact methods will become practical. 

We recently applied the truncated Wigner method to the 
problem of colliding condensates [7], showing that it was fea- 
sible to simulate a realistic three-dimensional system and pro- 
duce quantitative results. The aim of this paper is to expand 
upon and extend the treatment presented in there. To this end 
we provide here: 

a) A detailed derivation of the truncated Wigner method 
for colliding condensates, including a heuristic demon- 
stration of how we can justify the neglect of terms aris- 
ing in a Wigner function treatment of the problem, set- 
ting up the foundation of the method. We show in 
Sects. II B 2-IIB 3 that for this kind of problem, the va- 
lidity criterion for the method is that the density of the 
condensate (in co-ordinate space) must be very large 
compared with the density of the added quantum fluc- 
tuations. 

b) A treatment of the computational aspects of the prob- 
lem, which are rather subtle. The main feature to note 
is that the cutoff necessary in the effective field theory 
cannot be simply provided by the fineness of the spatial 
grid used for computations. Rather, in order to avoid 
aliasing, it must be provided explicitly by means of a 
projector. 

c) The evaluation of averages using full ensemble compu- 
tations. In [7] we evaluated averages using single com- 
putational runs, and averaging over regions of space 
where symmetry indicated the physics was the same. 
The results of our ensemble methods can provide addi- 
tional information on coherence properties of the final 
states. 

d) Results which can be compared with feasible exper- 
iments. At present there are no experimental data 
which can be quantitatively compared with the results 
of this paper. The work was motivated by the obser- 
vation of a strong halo in [3], but no detailed measure- 
ments were made on this halo which could be compared 
with theory — furthermore, the parameter regime is not 
fully within the range of validity of our methodology. 
The work reported in [24] has adapted our theoretical 
methodology to analyze a related experiment, and got 
good agreement. However, the calculation done was 
only two-dimensional, and therefore can only be re- 
garded as indicative. 

For this paper we have chosen a parameter regime 
which is experimentally attainable, which is in the re- 
gion where strong Bosonic stimulation is important, 
and which is fully within the region of validity of our 
methodology. In Sect. IV we give the values of appro- 
priate parameters for both sodium and rubidium con- 
densates. The numerical results are presented in SI units 



for a sodium condensate, and should be directly verifi- 
able experimentally. 



II. TRUNCATED WIGNER METHOD 

In dilute Bose gases the appropriate Schrodinger picture 
Hamiltonian is 
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H =j dx¥ (x) j- 

+ \J dx 1 ¥ (x') U 2h (x - x') <J (x')| * (x) . (1) 

Here the external potential is C/ cxt (x) and pairwise inter- 
actions between the bosons are characterized by the two- 
body scattering potential Uib (x — x'). The second-quantized 
field operator ^ (x) annihilates a particle from position x 
and obeys the equal time commutation relations for identical 
bosons 



*(x),#(x') = *t( x ) s $t( x ') 







*(x),* t (x') = <5(x-x'), (2) 
where 5 (x) is the three-dimensional Dirac delta function. 

A. Effective field theory 

1. Restricted basis field 

We now decompose the field operator onto a single-particle 
basis 
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where the mode operators obey the usual bosonic commuta- 
tion relations 
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By choosing the basis set to be the orthonormal eigenstates of 
the non-interacting portion of the Hamiltonian, Eq. (1), i.e. 
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the full Hamiltonian, Eq. (1) can be rewritten as 



H = Y^ fiujataj + ^Y1 ti r \ u 2b\ st) a]a, 
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It is usual to simplify the Hamiltonian, Eq. (6), by using an 
effective field theory, obtained by eliminating higher energy 
modes, whose time-dependence is so rapid as to be unobserv- 
able in experiments on ultra cold gases. This kind of proce- 
dure has a long history, and takes many different forms. The 
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description most appropriate to our methodology was given by 
Morgan [27], and divides the modes j into two sets, low- (L) 
and high-energy (H ) subspaces depending on whether hujj is 
less than or greater than a certain boundary energy e cut . Pro- 
vided e cut is sufficiently small, the effective Hamiltonian de- 
scribing the M low-energy modes can be found from Eq. (6) 
to be 



fl e ff = f^jOjClj 



U 
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where the interaction parameter is defined as U$ = AnfPa/m, 
for the s-wave scattering length a. 



2. Validity of the effective field theory 

In order for fl c ff to accurately describe the low-energy sub- 
space, e cut must be chosen so that the evolution of the high- 
energy modes is rapid compared to the evolution of the low- 
energy system, and L must include enough modes to ade- 
quately describe the system dynamics. For the colliding sys- 
tems we treat here, the low-energy subspace must be suffi- 
ciently large to contain all the modes required to represent the 
two colliding condensates and the scattering halo, as well as 
all the modes that can be directly involved in scattering events 
with those (condensate and halo) modes. In principle, this 
may conflict with the requirement that e cut be small — how- 
ever in the cases we consider here, the error is less than a few 
percent, as we describe in section IV. 



3. Projector representation 

It is useful to have some formal manner of decomposing ar- 
bitrary objects, such as field operators or wavefunctions, into 
components on either side of the boundary energy. To this end 
we define the orthogonal projection operators (projectors) 
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where V and Q act, respectively, as projectors onto the low- 
(L) and high-energy (H) subspaces. As these two subspaces 
completely span the infinite mode space, the projectors V and 
Q satisfy the closure relation V + Q = 1. 

In its coordinate space form, the low-energy projector acts 
on the arbitrary function / (x) as 

v [/ (x)] = ]T ^ (x) f dxV* (x') / (x') . (9) 

oeL J 

Applying this to the field operator given by Eq. (3) returns the 
restricted field operator 
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which is the component of the total field operator acting 
within the low-energy subspace. Because of the restricted na- 
ture of fy-p, the commutation relations given by Eq. (2) no 
longer apply. Rather it can be shown, using the mode opera- 
tor commutation relations (4), that the restricted field operator 
obeys the equal time relations 
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(x) , ¥ v (x')] = Mx,x'), 
where we have defined the restricted delta function 
Mx,x') = 5>*(x')^-(x). 
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Unlike the true (Dirac) delta function, the restricted delta 
function is spatially nonlocal, where the range of this non- 
locality scales as e c ^/ 2 . 



B. Wigner function evolution 

Let us define the density operator of the restricted basis field 
to be p (t), whose time evolution, using Eq. (7), is given by 
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The formulation of the truncated Wigner method is made 
using a multimode Wigner function representation of the den- 
sity operator p (t). The full details of how this is used are 
given in [26], and in brief are as follows. For a single mode, 
the Wigner function W (a, a* , t) is defined in terms of the 
Wigner characteristic function 

Xw (A,A*) = Tr{pcxp(Aa t -A*a)}, (15) 

as a Fourier transform 

W(a,a*) = -^ J d 2 \ exp(-Aa* + A*a) xw(A, A*). 

(16) 

The Wigner function exists for any density operator, and its 
moments are equal to those of the symmetrized operator prod- 
ucts: 



sym 



d 2 aa r (a*) s W (a,a*) (17) 
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The Wigner function is not guaranteed to be positive, but often 
is, and in these cases it behaves like a probability distribution 
for the variables a, a* . 
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The action of the mode operators a, a< on the density op- The extension to many modes is straightforward, 
erator can be expressed as the action of differential operators 
on the Wigner function using the operator correspondences, 
which follow from Eq. (16) 

ap{t) <-> (^a+^^JW(a,a*,t) (19a) 

l_c 
2 9a 
L_9 
2da 
l_o 
2 da 



P {t)a) ( a * + l--^]w{a,a*,t). (19d) 



For the restricted basis von Neumann equation, Eq. (14), we find using these operator correspondences that the evolution of 
the multimode Wigner function W ( { a.j , a* } , t) is given by 
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This Wigner function evolution is exactly equivalent to the 
von Neumann equation, Eq. (14). 



1. Wigner truncation 

Equations of motion for the Wigner function of the form 
given by Eq. (20) are well known, particularly in quantum op- 
tics [26], and even in the case of a few variables are not easy 
to solve numerically. The Wigner truncation, in which the 
third-order derivative terms are dropped, has often been made, 
since the resulting equation of motion, having only first-order 
derivatives on the right, is of the form of a Liouville equation. 
It thus describes an ensemble of trajectories obeying an equa- 
tion of motion which is essentially classical, and which can be 
simulated relatively straightforwardly. 

The justification for this truncation is intuitively reasonable; 
if the quantum state of the system is such that it is "almost 
classical", then the classical equation should prevail. One 
can present scaling arguments, in which it is assumed that 
the Wigner function behaves like a sharply peaked probabil- 
ity distribution centered around a macroscopic mean value of 
the parameters a,j . These arguments can then be used to show 
that the contribution from the third-order derivative terms is 
negligible provided the mean values of all of the aj are large. 
This kind of argument has been made relatively rigorously by 
Polkovnikov [28], who has shown how the third-order deriva- 
tive terms give a correction to the classical trajectories. 

In the case of colliding condensates this criterion is not 
valid. The initial state of the system is only highly occupied 



in the modes in the vicinity of the two incoming momenta, 
and we want to consider the evolution into a large number 
of initially unoccupied modes. However, experience in quan- 
tum optics shows that the Wigner truncation can be valid for 
such a system; the prime example is the degenerate paramet- 
ric oscillator, in which two modes of the electromagnetic field 
are made to interact by means of a nonlinear crystal to give a 
Hamiltonian of the form 

H dpo = huja { a + 2huPb + g |b (a f ) 2 + S f a 2 | . (21) 

Here a and b are destruction operators for field modes with 
frequencies u and 2u. A Wigner function treatment soon re- 
veals an equation of motion with first- and third-order deriva- 
tives with respect to corresponding Wigner function variables 

a, a* , (3, (3* . 

In the relevant physical situation, the mode a is initially 
unpopulated and the mode b is derived from an intense laser, 
and as such is essentially a classical field — thus one makes 
the replacement 

b^Be~ 2iu] \ (22) 
where B is a classical amplitude, and thus 

ffdpo -> ffdpo = f^a)a 

+ g^Be- 2lult (ajf + B*e 2lult a 2 Y (23) 

This approximate Hamiltonian produces a Wigner function 
equation of motion with only first-order derivatives with re- 
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spect to a, a*, and is well verified to give accurate physical 
predictions for the production of quanta in the mode a. 

Alternatively, one could simply drop the third-order terms 
in the Wigner function equation of motion, and one would get 
equivalent results in the limit that one could neglect depletion 
of the field b. This is in spite of the fact that the occupation of 
the mode a is not large — thus it appears that it is sufficient 
for only some of the modes to be highly occupied for it to be 
valid to drop the third-order derivative terms. 



2. Validity of the Wigner truncation for colliding condensates 

The case of colliding condensates is analogous to that just 
discussed. There are a number of highly occupied modes — 
those corresponding to the original condensate packets, and a 
larger number of unoccupied modes. In the following we give 
a heuristic analysis of why the Wigner truncation should be 
acceptable for this system. 



Consider a multimode Wigner function, which at some time 
t is factorizable into single-mode functions, of the form 



W ({ay, a}}, r) = [] — ex P 
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where ctj gives the expectation value (coherent) amplitude of 
the jth mode, and Tj < 2 is the inverse width of the Wigner 
function. This Wigner function includes both thermal and co- 
herent (for which Tj = 2) statistics, but does not allow the 
modes to exist as pure number states or other more elaborate 
forms. In fact we use exactly this Wigner function when con- 
structing the initial states of our simulations, for which we 
find that the modes display essentially Gaussian statistics at 
all times, with minimal correlations between the modes. Thus 
we expect the following analysis to be appropriate for the du- 
ration of the collisions considered here. 

Substituting the Wigner function given by Eq. (24) into the 
nonlinear portion of Eq. (20) gives the evolution at time r as 
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Here the first bracketed set of terms within the braces arise 
from the first-order derivatives, while the second bracketed 
set of terms arise from the third-order derivatives. 

Analogously to the definition of the restricted basis field 
operator, as given by Eq. (10), we define the classical wave- 
function 



*p (x) = $i ( x ) 1 



'J- 



(26) 



which represents a possible state of the total restricted basis 
field, including both the condensate and noncondensate parti- 
cles, at any given time. We also define the related wavefunc- 
tion 



jet 1 



(27) 



whose expectation value, calculated using the Wigner func- 
tion given by Eq. (24), is found to be 



fr (x) = (x)}, 



EtWj«i.- (28) 
jeL 1 



Using these wavefunction forms in the Wigner function evo- 
lution, Eq. (25), gives 



BW ( nonlin J r 
iK-ft = 2C/o J dx [ (Cv - Cv ) *p - n - fa )] 
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where we have suppressed the explicit spatial dependences and have retained the ordering of Eq. (25). 
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To justify the Wigner truncation, we now show that the 
terms arising from the cubic derivatives in Eq. (29) are small 
compared with the first order derivative terms for all points x 
on the coordinate space field. This local analysis requires that 
the inequality 



1*7 



« 1, 



(30) 



should hold over all space. A useful quantitative description of 
the inequality is in terms of the distributional averages (clas- 
sical expectation values), for which we find that 



I6p - e 
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which is of similar magnitude to 5p (x, x), Eq. (12). For the 
expectation value of the \^f-p\ 2 however, a more useful repre- 
sentation can be obtained using the correspondence of Wigner 
function averages to the quantum expectation values, Eq. (17). 
We find that in the general case (i.e. irrespective of the partic- 
ular form of the Wigner function) 



(x)| 2 )^=n(x) + i^(x,x), 



(32) 



where the total density of real particles n (x) is defined using 
the restricted basis field operators by 



,(x) = (n(x)) (x)%, (x) 



(33) 



Using Eqs. (31,32), the validity criterion for the Wigner trun- 
cation, Eq. (30), becomes 



For a zero-temperature homogeneous field, such that Tj = 
2 for all modes, the validity condition becomes simply N ^> 
M, which is similar to that given by Sinatra et al. [29]. How- 
ever, for the inhomogeneous finite-temperature case, the lo- 
calized truncation condition given by Eq. (34) is less easy to 
justify, especially when one considers that for these inhomo- 
geneous fields there may be regions where the total particle 
density goes to zero. However, the part of the Wigner function 
evolution dependent upon interparticle scattering is significant 
in only those regions where there is a high particle density, i.e. 
the regions where the truncation is accurate, so that the trun- 
cation can be made over all space without adversely affecting 
the accuracy of the approximation. Interestingly, this justi- 
fication relies on the relative magnitude of the particle and 
mode-function densities, rather than the numbers of particles 
and modef unctions. Thus it appears to be possible to accu- 
rately apply the Wigner truncation to systems in which there 
are significantly more basis modes than real particles. 

3. Critique of the Validity Criterion 

The validity criterion in the form Eq. (34), or in the form 
for a homogeneous system N ^> M, shows that for a given 
number of real particles in the system, accurate results will 
not result if the number of basis modes M is too large. This 
is fundamental to the truncated Wigner function method, in 
which vacuum noise is added to every mode. Methods based 
on the P-function[l 1, 12, 26] do not have this problem, since 
noise is added dynamically and only to modes with real occu- 
pation. However, as noted in the introduction, other technical 
difficulties so far make these methods more difficult to use in 
practice. 
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Thus in order to justify the truncation it is required that the 
real particle density is large compared with the commutator 
of the restricted field, dp (x, x). 



4. Truncated Wigner function Fokker-Planck equation 

Within the validity regime of the Wigner truncation the 
Wigner function evolution given by Eq. (20) is well approxi- 
mated by the truncated Wigner function Fokker-Planck equa- 
tion 



ih 



dW 
~dt 



E I twjOtj + t/o ^ (jr\st) a*a s a t 1 W 
jeL i I rsteL ) 

- E ^ I + u ° E 0' r i st > 
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a r a* s al > W. 



(35) 



Note that we have also removed from Eq. (35) the nonlin- 
ear evolution dependent upon 5 r . s (see Eq. (20)), which cor- 
responds to terms dependent upon Sp (x, x) in coordinate 
space, and whose influence on the total evolution, as we dis- 
cussed above, is negligible. Indeed, in order to preserve en- 



ergy conservation these terms must be removed alongside the 
cubic derivative terms. 
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C. Ensemble differential equations 

The Liouville equation, Eq. (35), gives the equation of mo- 
tion for the ensemble of trajectories of the variables { otj , a* } . 
The corresponding equations of motion for individual trajec- 
tories of the system are given by [30] 



ih—r- 

dt 



a*a s a t , (36) 



rsteL 



where j e L. Every distinct realization of the differential 
equation uses a different initial noise field, and hence yields 
a distinct trajectory in the phase space of the system — we 
describe appropriate initial states in section II C 2. 

By using our previously defined restricted basis wavefunc- 
tion, Eq. (26), in the time-dependent form ^-p (x, t), we can 
rewrite the evolution of the jth low-energy mode amplitude as 



dt 3 3 



+ U / dx^l^r*^, (37) 



which we find to be convenient for numerical simulation. We 
use Eq. (37) (in a dimensionless form) to calculate the central 
results of this paper. 

Note that although our differential equations are not 
stochastic differential equations, they do have a stochastic na- 
ture that arises from the random component of the initial field. 



1. Projected form of the differential equations 

It is instructive to obtain the differential equation for the en- 
tire restricted field, rather than the evolutions of the individual 
mode amplitudes. Again using the definition of the restricted 
basis wavefunction, we find that Eq. (37) gives rise to 



iH 
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(38) 

where we have replaced the basis eigenenergies with the cor- 
responding diagonal Hamiltonian and have recognized the 
form of the projector onto the low-energy mode space, Eq. (9). 

It is straightforwardly shown from any of Eqs. (36,37,38) 
that the normalization of the total field for a single trajectory, 
defined as 



Af = 



| 2 = /dx|* 7 ' 2 



(39) 



is strictly conserved, so that dN /dt = 0, as is the total field 
energy 



JdL 



Uo 
2 



dx l^-pl 



(40) 



Note that JV and £ are not the physically observable field pop- 
ulation and energy, as they include contributions from the vir- 
tual particles (the noise). It could be imagined that the pro- 
jector serves to remove those quanta that are shifted outside 



the low-energy subspace by pairwise collisions, reducing the 
total population of the field over time. This view is incorrect. 
Rather the projector simply disallows these processes, a result 
that is most easily seen from Eq. (36). 



2. Initial states 

As the ensemble differential equations contain no dynamic 
noise sources, to sample the evolution of the Wigner func- 
tion we are obliged only to ensure that the ensemble of ini- 
tial states yields the appropriate Wigner function. For the 
zero-temperature collisions presented here we assume an ini- 
tial Wigner function of the form given by Eq. (24), where the 
modes are assumed to display uniformly coherent statistics, 
for which Tj = 2 for all j. The corresponding initial state of 
the jth mode amplitude for a single trajectory is given by 



aj (0) 




(41) 



Here, as before, ctj is the coherent amplitude of the mode, 
and can therefore be identified as the condensate amplitude 
for the jth mode. The quantities Aj,Bj are real independent 
Gaussian random variables of zero mean and unit variance, 
such that 



(A,) = (Bj) = {MB,) = 



(AiAj) = (BiBj) = 5. 
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(42) 



The equivalent coordinate space initial state is found using 
Eqs. (26,41) to be 



*7> (x,0) 



Y^Tpj CXjo (0) 

^ Vo (x, 0) + xv (x) 




(43) 



(44) 



Here ^-p = (^-p)w ' s me coneren t field amplitude, i.e. the 
condensate wavefunction. The spatial fluctuations, xv are 
found from Eq. (42) to satisfy 



(Xv) = 0, 



(Xv (x')xp (x)) = 



■</>*(*') <Mx), 



(45) 
(46) 



where Eq. (46) can be seen to be similar to the restricted delta 
function defined by Eq. (12), and proportional to it when all 
Tj = 2. For our uniformly coherent initial state we find that 
the expected total field normalization, as given by Eq. (39), is 



(AO=iV (0) + ^, 



(47) 



where 



N (t) = K W\ 2 = I dx |*t>o (x, t)\ 2 , (48) 

3&L J 

is the total number of condensate particles and, as before, M 
is the number of low-energy modes. Although the number 
of condensate particles is in general not conserved, the total 
number of real particles N is conserved. 
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III. NUMERICAL METHODS 
A. Computational units 

It is convenient for the purposes of numerical simulation to 
express the equations in dimensionless computational units. 
The systems we treat are initially confined within a harmonic 
potential 



E/harm (x) = — {u^X 1 + UJ 2 y y 2 + W 2 Z 2 ) , 



which provides 



to — — , 



(49) 



(50) 



as a natural choice of units for length, time and energy respec- 
tively. 

In these units Eq. (37), becomes 



.dotj 
1 di 











u> jaj - 


hft j 


f dx$ 





(51) 



where we have used tildes to indicate quantities in compu- 
tational units (ctj is identical in both computational and S.I. 
units). As an example, the dimensionless interaction parame- 
ter U a is 



U 



Up 



8™ 



2mu! r 



(52) 



B. Plane wave basis 

While our formalism allows the use of any set of orthonor- 
mal basis states tpj (x), the most appropriate choice of basis 
for collisions occurring in free space is the plane wave states. 
We use a three-dimensional, periodically-bounded rectangular 
space of volume V = L x x L y x L z , for which the (normalized 
to unity) plane-wave modes are 



^■(x) = -l= e ^- x . 
Here the wavevector k_, associated with the jth mode is 



L x 



2imj 



L z z 



(53) 



(54) 



where mj, rij andpj are integers. 

In our dimensionless computational units, the single- 
particle energy of the jth mode is given by uJj = fc 2 , and the 
differential equation describing the evolution of the jth mode, 
Eq. (51), becomes 



daj 
i — 
dt 



~k)a 3 + 



U 



dice 



(55) 



We have previously defined the low-energy mode subspace 
L to consist of all those modes whose energies are less than 



some cutoff energy e cut . It is easily seen that for the plane- 
wave modes this energy cutoff is translated to a spherical cut- 
off in wavevector space. Thus the L subspace consists of all 
those modes whose wavevectors satisfy kj < k cut , where 
the cutoff wavenumber is defined (in computational units) as 

^-cut = V^cut' 



C. Propagation algorithm 

To propagate the differential equations we employ a modi- 
fied version of the Fourth-order Runge-Kutta in the Interac- 
tion Picture (RK4IP) algorithm [31], which has been used 
extensively to simulate the time-dependent Gross-Pitaevskii 
equation. The only difference in our equations of motion is 
the presence of the projector, which is straightforwardly per- 
formed in momentum space, albeit with a small number of 
required conditions that are worth stating explicitly. 

Using a plane-wave basis, the conversions between coordi- 
nate space and mode space required by Eq. (55) are achieved 
using Fast Fourier Transform (FFT) algorithms [32], which 
require that the two spaces are represented using identical 
numbers of rectangularly arranged grid points. Thus, given 
that the low-energy mode space is spherically bounded, we 
must include additional modes (from the H subspace) to form 
a rectangularly bounded space (the padded mode-space). Us- 
ing a padded mode space that tangentially bounds L is not 
adequate, because with such a grid population scattered from 
the L subspace into the H subspace, which should be removed 
by the projector, can be aliased by the FFT back into the L 
subspace, thereby adversely affecting the accuracy of the sim- 
ulation. A rectangularly bounded mode-space grid of extent 
4A: C ut x 4A: cu t x 4fc cut is the smallest grid that prevents aliased 
components returning to the L subspace, and therefore allows 
for accurate calculation of the low-energy subspace dynamics 
while minimizing computational memory requirements. This 
padded mode-space has ~ 48 /it times as many grid points 
as there are low-energy modes, and significant computational 
memory savings can be made by representing the L subspace 
on the full padded grid only when absolutely necessary for the 
propagation algorithm. 



IV. SYSTEM OF INTEREST 

In our previous letter [7] we presented results obtained from 
a single trajectory describing a very similar collision to one 
that had been realized experimentally [33]. However, care- 
ful analysis indicates that the systems presented in [7] may 
stretch the validity criteria of the truncated Wigner method, 
in particular due to an inappropriately low k cu t- In this paper 
we treat colliding systems that probe a different region of pa- 
rameter space, one which is both experimentally possible and 
well within the validity regime of our formalism. 

The condensate part of our initial state (see Eq. (44)) is as- 
sumed to be composed of two equally populated wavepack- 
ets derived from a single harmonically trapped Bose-Einstein 
condensate using a short it/2 Bragg pulse [34] of wavevector 
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(Ag, 0,0), so that 



0) = ^(x) 



(56) 



where the envelope function (x) is that of the initial, un- 
scattered condensate wavefunction. We describe this enve- 
lope wavefunction by the No atom ground state solution to 
the time-independent Gross-Pitaevskii equation. Eq. (56) as- 
sumes a centre-of-mass frame, which reduces the number of 
low-energy modes required and provides a convenient sym- 
metry. We assume a zero-temperature initial state, such that 
the all modes are initially coherent and N = No (t — 0), and 
we remove the confining potential at t = 0. 

We use a (dimensionless) nonlinear parameter of Uo = 1 x 
10 -2 , an initial (total) condensate population of two million, 
and a relative wavepacket wavenumber of Aq = 10. This 
parameter set corresponds to a system of 23 Na atoms, initially 
confined within a potential with u> x = 2tt x 4.57 Hz, where 
the wavepackets move with a relative speed of 4.0 mms -1 , 
or alternatively to a 87 Rb system initially confined within a 
trap with ui x = 2ir x 0.31 Hz, where the wavepackets are 
separating at 0.53 mms -1 . We assume the trapping potential 
is cylindrically symmetric, with A = X z /X x , y = a/8. 

For our (dimensionless) parameter set, the Thomas-Fermi 
chemical potential [35] is found to be /ixF = 21.4, giving 
Thomas-Fermi radii of xtf = Vtf = 2\/fi-TF = 9.24 and 
ztf = 2a//2tf /A = 3.27. In order to enclose the colliding 
system until packet separation we find that a volume of V = 
48.9 x 33.5 x 33.5 is appropriate. 



A. Truncated Wigner method validity criteria 

In restricting the system to the low-energy subspace, as 
described in section II A 2, it was stated that the momentum 
space must be large enough to include all relevant modes 
to describe the evolution. For the colliding system this re- 
quirement translates to k cut > 3Aq/2, such that all possi- 
ble scattering events directly involving the initial condensate 
wavepackets and the halo are included. To meet this validity 
criterion we use a low-energy subspace cutoff of A,- cut = 18, 
for which M = 5.4 x 10 6 . Using a result from [36] we find 
that the error in the s-wave scattering length in the effective 
Hamiltonian, Eq. (7), is approximately 0.5% for this value of 
the cutoff. 

The second validity criterion for the truncated Wigner 
method, Eq. (34), is most easily expressed as 



n (x) » Sp (x, x) 



(57) 



For a harmonically trapped Thomas-Fermi condensate, the 
(dimensionless) maximum density is h (x = 0) = fiTv/Uo, 
so that for our system n (x = 0) = 2140. For a plane wave 
basis we find that the equiposition restricted delta function can 
be expressed as 



S-p (x,x) 



M 

V 



p 

6tt 2 ' 



(58) 



For fc cut = 18, we find that S-p (x, x) ss 100, so that we expect 
the Wigner truncation criterion, Eq. (34), to be satisfied for 
this system. 

We have simulated (otherwise identical) additional trajec- 
tories using a larger cutoff, fc cut = 20.57, for which M = 
8.4 x 10 6 . Using the extrapolation method outlined in section 
V B 3, we find that the difference in the calculated total co- 
herent population between the system with this larger cutoff 
and our principal system approaches 1.25% by the end of the 
simulation times. 

For the remainder of this paper we present all results in SI 
units, appropriate for the 23 Na system described above. 



V. RESULTS 
A. Single trajectory results 

1. Momentum space 

In Fig. 1 we show the mode populations for a single trajec- 
tory of our colliding system, for those modes whose velocities 
lie on the planes v z = and v x = 0, at a sequence of times. 
There are two major points of interest in this figure. First, 
we observe from the v z = planes that the original conden- 
sate wavepackets broaden and change shape from ellipsoidal 
to crescent shaped. Secondly, we observe the generation of a 
circular feature, centered at the system centre of mass veloc- 
ity; this is the scattering halo that is the focus of this work. 
Note that the mode population scale of Fig. 1 is logarithmic, 
where the scale has been chosen to display the lower popula- 
tions to best effect. A consequence of this is that the higher 
populations saturate for \ctj\ > 100, which is rather lower 
than the population of 3.2 x 10 4 for a mode at the centre of 
each of the condensate packets at t = 0. 

We can see from Fig. 1 that the scattering halo first appears 
(in our centre of mass frame) centered at a radius slightly less 
than the initial central wavenumber of the individual conden- 
sate wavepackets, due to the interaction energy cost associated 
with creating perturbations in the field [5]. As time progresses 
the halo broadens, largely inwards, to occupy (in the average) 
virtually all those modes within a certain wavenumber radius. 
This broadening is a consequence of scattering events between 
particles already present in the halo with those in one of the 
condensate wavepackets, or other halo particles. Both pro- 
cesses serve to redistribute population within the halo. In ad- 
dition we note from Fig. 1 (c) that when the halo first appears 
there is an angular dependence, with those modes whose po- 
lar angles relative to the collision (v x ) axis are closest to ir/2 
having the greatest increase in population. 

This anisotropy occurs despite the inherent isotropy of s- 
wave scattering, because rather than the scattering of single 
quantum, here we are dealing with the collective scattering of 
many particles in the presence of a matter-wave grating. The 
probability of any single scattering event is thus proportional 
to the overlap integral between the input and output wave- 
functions, which means that modes whose momenta are either 
perpendicular or parallel to the direction of the grating have 



10 



(a) 


(b) 




(d) 


• • 




(e) 


if) 

-' 


r i i 






(h) 

■» 


-2 2 

[lllffiS -1 ] 


-2 M 2 
*?a [mms -1 ] 



10 



-U.5 i n*'.l 



10 J.r. 10 2.ii 



in 



(a) 


(b) 


■I | | I 


w t 


-4 4 


A 4 


IV [mms -1 ] 


v x [mms -1 ] 




-o.b 10 o.o m m 


;o 10 taM io 



FIG. 2: Single trajectory velocity mode populations on the plane 
v y — at 1.6 ms (a) and 37.7 ms (b) for the same collision as in 
Fig. 1. The low-energy subspace boundary is visible as a circle, be- 
yond which no population occurs. 
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FIG. 3: Single trajectory velocity mode phases on the plane v z — 
at 8.2 ms (left) and 16.4 ms (right), corresponding to Fig. 1 (c,e). 



FIG. 1: Velocity mode populations on the planes v z = (a,c,e,g) 
and v x = (b,d,f,h) at times t = (a,b), 8.2 ms (c,d), 16.4 ms (e,f) 
and 37.7 ms (g,h) for our colliding system described in the text. Note 
that the simulation volume exceeds that shown. 



the largest overlap integral and hence experience the greatest 
growth. 

Although not visible in Fig. 1, three- wave mixing between 
the condensate packets gives rise to additional wavepackets 
centered (in momentum space) at k x = ±3Ag/2. However, 
due to the kinetic energy mismatch involved in their genera- 
tion, these packets are transient, and have essentially disap- 
peared by the time the initial wavepackets are separated. This 
effect is shown in Fig. 2, where these higher order packets 
are present shortly after the collision begins, and are absent 
by 37.7 ms. Additionally, this figure shows the asymmetry in 
the condensate wavepackets both initially and as a result of 



anisotropic broadening due to the oblate nature of the initial 
trapped condensate. 

The scattering halo is not uniformly populated. Rather it is 
characterized by distinct patches of high population separated 
by regions of low population. In Fig. 3 we plot the phases of 
the modes on the plane v z = at times corresponding to the 
second and third rows of Fig. 1. From Fig. 3 (left) we observe 
that as the population in the halo begins to grow, small re- 
gions of relatively constant phase are formed, which by com- 
paring Figs. 1 and 3 can be identified with those regions of rel- 
atively high population. Hence we label these small, highly- 
populated regions of momentum space phase grains. The size 
and aspect ratio of these phase grains is very similar to those 
of the initial condensate wavepackets in momentum space, a 
result that we explore more quantitatively in section V E. We 
note that the phase profile of each grain is established prior to 
any significant population gain. 
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m 




(0 




tinct types of wavefunction, a relatively smooth outer portion 
and a fragmented central region. In section VB, we iden- 
tify the smooth outer shells as condensate that survives the 
collision, while the central region forms the scattered halo, 
together with a very small condensate remnant. The high- 
est density during the course of the simulation remains close 
to the origin, which somewhat counterintuitively leads to the 
scattered quanta being more dense than the condensates at 
later times. After separation of the condensate packets, scat- 
tering events taking place within the central (fragmented) re- 
gion involve at most one condensate quantum, and redistribute 
population within the halo. 

Estimating the separation time using Thomas-Fermi con- 
densate wavepackets of unchanging shape returns 30.6 ms, 
significantly longer than that observed in our simulation. The 
reduction in separation time reflects the significantly distorted 
condensate packets, from initially ellipsoidal to a shape most 
easily described as a hemi-ellipsoidal shell (see Fig. 4 (g)). 

We observe from Fig. 4 (g) that near the end of the simu- 
lation, the coordinate space field is beginning to resemble the 
momentum space field, as shown by Fig. 1 (g). By this time 
the rate of change of the mode populations is very small, and 
we therefore expect that Fig. 1 (g) will give an excellent in- 
dication of the coordinate space density distribution for the 
system following a sufficiently long ballistic expansion. 
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FIG. 4: Coordinate space density for the same collision as Fig. 1 on 
the planes z — (a,c,e,g) and x = (b,d,f,h) at times t = (a,b), 
12.6 ms (c,d), 25.1 ms (e,f) and 37.7 ms (g,h). Note the logarithmic 
density scale. 



2. Coordinate space 

In Fig. 4 we show the evolving total coordinate space den- 
sity |\T/-p (x, t)\ 2 on the planes z — and x — for the tra- 
jectory shown in the preceding section. The initial state of the 
field, shown by Fig. 4 (a,b), displays the modulated ground 
state condensate profile given by Eq. (56), while the vacuum 
noise component is observed to be spatially uniform in its av- 
erage density, as required for the plane-wave basis, extending 
into and distorting the condensate profile. 

As the collision proceeds we observe that at the centre of 
the system, where the particle density is highest, the regular 
fringe pattern arising from the overlapping wavepackets be- 
gins to break down, indicating that quanta are being generated 
with momenta other than those of the condensate wavepack- 
ets. By 25.1 ms into the collision, as shown in Fig. 4 (e), the 
fringes appear to be completely absent. Following this break- 
down, the high density region of the field comprises two dis- 



3. Turbulence 

Within the scattering halo, a large number of vortices have 
been detected between the highly populated grains in both ve- 
locity and coordinate spaces. We demonstrate this in Fig. 5, 
where we show the velocity mode populations on the plane 
v z = and the coordinate space density on the plane 2 = 
at the end of the collision, together with the detected vortices. 
Note that we plot in Fig. 5 only those vortices that reside in 
regions of relatively high average mode population (Afj > 1) 
or density (l^-p | 2 > 3 x 10 11 cm -3 ), which filters out those 
vortices that are present in the field solely by the presence of 
the virtual particle fluctuations, and those that are physically 
observable, although the threshold choice is somewhat arbi- 
trary. 

We have observed (but not plotted) that a significant growth 
in the number of physical vortices coincides with the growth 
of the scattering halo, indicating that, in much the same way 
as the halo growth can be viewed as amplification of the vac- 
uum fluctuations, the observed vortices are amplifications of 
the underlying quantum vortices. This provides the basis for 
our identification of quantum turbulence within the scattering 
halo. 



B. Coherent and incoherent fields 

We have assembled an ensemble of 80 individual trajecto- 
ries of our colliding system, where the initial states of each 
differ only in the applied virtual particle noise, as described in 
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FIG. 5: (a) Velocity mode populations on the plane v z = and 
(b) coordinate space density on the plane z = at 37.7 ms into 
the collision. Crosses and circles respectively show unit vortices of 
positive and negative sense. 



section II C 2. In the following sections we use this ensemble 
to calculate various quantum statistics of the colliding system. 



1. Momentum space 

Using the correspondence between moments of the Wigner 
function and symmetrized operator products, Eq. (17), we can 
calculate the total expectation population of the jth mode, 
Nj (t), using 



N j (t) = (N j )(t) = (a]a j )(t) = (\a j (t)\ 2 ) 
Using the coherent amplitude of the jth mode 

a i0 (*) = (%) 0) = («j (*)>W! 



1 

2' 
(59) 

(60) 



we can calculate the coherent (i.e. condensate) part of the total 
mode population as 



Njo (*) = |(«j(*)> 



irl 



(61) 



In Fig. 6 (a,c,e,g) we plot the coherent mode populations 
on the plane v z = for our colliding system. From these 
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FIG. 6: Coherent (a,c,e,g) and incoherent (b,d,f,h) velocity mode 
populations on the plane v z = at times (a,b), 8.2 ms (c,d), 
16.4 ms (e,f) and 37.7 ms (g,h). 



plots we observe that the condensate population is restricted 
to two wavepackets that, consistent with the behavior shown 
in Fig. 1, broaden and change shape over the course of the 
collision. A benefit of plotting the coherent populations only 
is that the small-scale structure of the condensate packets is 
rather more apparent here than in Fig. 1. Although not shown 
in Fig. 6, the higher order wavepackets, when they appear, are 
also found to contain condensate particles. 

The corresponding incoherent mode populations are calcu- 
lated as the difference 



N 



(incoh) /,n _ 



(t) = Nj (t) - N JQ (t) . 



(62) 



In Fig. 6 (b,d,f,h) we plot the incoherent mode populations 
on the plane v z = for the collision. From these plots we 
observe that from an initial state with zero incoherent popu- 
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lation (as required by our assumed initial Wigner function), 
incoherent population builds up initially as a semi-spherical 
scattering halo, where those modes occupied by the conden- 
sate wavepackets have a relatively small incoherent popula- 
tion. Subsequently population is transferred, via further scat- 
tering events, to broaden the distribution of incoherent quanta. 



2. Coordinate space 

The coherent and incoherent fields can also be calculated in 
coordinate space. Using the coherent probability amplitude 



(x,i) = (*7> (x))(i) = <*7> (X,i)) 



W ' 



(63) 



the coherent particle density (condensate density) can be cal- 
culated as 



n (x,i) = \V Vo (x,t)| 2 



( x ) a io (*) 



(64) 



The corresponding incoherent particle density is found using 
both the total particle density, Eqs. (32,33), and the coherent 
particle density to be 



(incoh) 



(x, t) = n (x, t) — n (x, t) 



(65) 



2\ 



*7>(x,t)| J 



W 



\(*v (x,t)> 



in 



2<M X , X ) 



(66) 



where the restricted delta function, S-p (x, x') is defined by 
Eq. (12). We plot both the coherent and incoherent particle 
densities for our colliding system in Fig. 7. 

From Fig. 7 (a,b) we observe that the initial state is uni- 
formly coherent, consistent with the assumed initial Wigner 
function. As the collision progresses we observe that signif- 
icant incoherent particle density develops close to the origin, 
generated from the condensate particles. Note that Fig. 7 (e,g) 
shows that a small interference pattern exists in this central re- 
gion even late into the collision, which was not apparent from 
a single trajectory (see Fig. 4). 



3. Total coherent and incoherent populations 

One of the important quantities that can be calculated from 
these simulated collisions is the total number of particles scat- 
tered out of the condensate wavepackets. In our previous Let- 
ter [7] we calculated this quantity approximately using a spa- 
tially dependent counting method in momentum space. Such 
methods are routinely used experimentally. However, using 
our ensemble we can obtain the total coherent population 



AT (t) = £K (t)r 



(67) 
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FIG. 7: Coherent (a,c,e,g) and incoherent (b,d,f,h) particle densities 
for our colliding system on the plane z = at times (a,b), 12.6 ms 
(c,d), 25.1 ms (e,f) and 37.7 ms (g,h). 



and incoherent population 

^(incoh) ( t ) = N- N a (t) , 



(68) 



with a great deal more accuracy. 

A well-known result for averages of randomly distributed 
variables is that the convergence of a finite number of samples 
M. to the true value scales as 1/ \J~M. (where convergence is 
measured as the error in the calculated average). We define 
the total coherent field population, calculated from Ai trajec- 
tories, as 



1 



M 



T7>>i(*) (m) 



M ^ 

m— 1 



(69) 



where otj (t)^ is the time-dependent mode amplitude of the 
?7ith trajectory. Given the large number of modes within our 
colliding system, and the magnitude squaring operation in 
Eq. (69), it is expected that 



N (t) 



(M) 



N (t) 



C(t) 
M ■ 



(70) 



where N (t) = N (t) (oo) . Here C (<) is some dimension- 
less time-dependent function, independent of M., whose ex- 
act form is in fact unimportant for our purposes. Assuming 
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FIG. 8: Total coherent (a) and incoherent (b) populations for the col- 
liding system calculated using the extrapolation method outlined in 
the text. 



that Eq. (70) holds, we can therefore use any two calculated 



N (t) 



(M) 



of dissimilar M. to obtain the ensemble limit. In 



particular, using M. = 1 for one of these two, we find that 



N (t) 



No(t) (1) ~MN (t) M 
1-M 



(71) 



In practise we find that this extrapolation method is extremely 
accurate using just two trajectories. This result makes a sur- 
vey of the parameter dependence of the global field quantities 
numerically efficient although, unfortunately, an equivalent 
extrapolation method for local field quantities is significantly 
less accurate. 

In Fig. 8 we plot the total coherent and incoherent popula- 
tions for our colliding system, calculated using our extrapola- 
tion method with M = 80. From these curves we observe that 
the coherent population monotonically decreases with time, 
with the most significant depletion occurring between 5 ms 
and 25 ms into the collision, by which time 30% of the total 
population is incoherent. Note that the incoherent population 
continues to increase even after separation of the condensate 
wavepackets, which results from decohering interactions in- 
volving individual wavepackets only. Note also that the dif- 
ference in the calculated populations is less than 1.5% for all 
M € [2,80]. 
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to the symmetrized quantum expectation values, Eq. (17), we 

.( 2 



find that gf (t) can be calculated using 



4 2) (*) - 



Wl 4 ) 



ir 



Moi(*)l a ) 



aj (t)\' 



(73) 



C. Local correlation functions 

1. Momentum space 

The normalized second-order equiposition mode-space cor- 
relation function 




probes the quantum statistics of the jth momentum mode. If 
a mode displays coherent statistics, then we should observe 

(2) 

that g - ' = 1. Conversely, if the mode displays Gaussian (i.e. 

(2) 

thermal) statistics, then we should find that = 2. By 
applying the correspondence of the Wigner function moments 



We plot gj (t) on the plane v z — for our colliding sys- 
tem in Fig. 9. Due to the normalized character of the corre- 
lation function and the finite number of ensemble members, 
those modes with small populations give highly noisy results. 
Therefore we have plotted in Fig. 9 only those modes whose 
real particle population is larger than one half. (Of course 
the quantum statistics of a mode are only defined when that 
mode is populated.) From Fig. 9 we observe that the conden- 

(2) 

sate wavepackets have g y - = 1, and are therefore coherently 

(2) 

populated, whereas the halo modes have <^ ; « 2, character- 
istic of Gaussian statistics. The higher level of noise in the 
results returned for the halo modes as opposed to the con- 
densate modes is a consequence of the much lower average 
population of the halo modes, and would be reduced for an 
increased number of ensemble members. 
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D. Pair creation correlation functions 

The halo formation process can be viewed as a four-wave 
mixing, in which the two input particles are from each of the 
condensate wavepackets, and the output particles appear (in 
the centre of mass frame) in modes of approximately (due to 
the finite momentum spread of the condensate wavepackets) 
opposite momentum. It is expected therefore that modes of 
opposite momenta will display correlations in both their am- 
plitude and population, at least at early times. 

An appropriate (normalized) correlation function to quan- 
tify the amplitude correlation between modes of opposing mo- 
menta is 



(0,2) 



(*) 



(a^) (t) ~ (t) (aj) (t) 



(76) 
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FIG. 10: Normalized second-order equiposition coordinate space 
correlation function g^ 2 ' (x, t) for the collision on the plane z — 
at times (a), 12.6 ms (b), 25.1 ms (c) and 37.7 ms (d). Results are 
only shown for those points where n (x) > 3 x 10 11 cm -3 . 



where we use k, 



-k; and the second term in the numera- 



tor corrects for the non-zero amplitude expectation value for 
condensate modes. This function can be written in terms of 
moments of the Wigner function as 



2. Coordinate space 

We can similarly characterize the quantum statistics of the 
system in coordinate space, for which the appropriate ana- 
logue to Eq. (72) is 



g {2) (x,t) 



v v (x) vip (x) m v (x) ( X ))(t) 



¥ v (x) i> v (x) ) (t) 



(74) 



By expanding the field operators on the restricted mode basis, 
Eq. (3) and applying the Wigner function moment correspon- 
dences, Eq. (17), we find that the second order normalized 
coordinate space correlation function can be calculated using 
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(75) 
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We plot the magnitude of g^j (t) for kj = — k^ at a se- 
quence of times in Fig. 11. (Due to our choice of i, j mo- 
menta the results are symmetric about the origin.) From these 
slices we observe that there is a definite amplitude correlation 
between modes of opposite momenta, largest at early times 

(5ij ' 2 ' > (t) ~ 0-3 at 6.6 ms), and decreasing as the collision 
progresses. The degradation in the correlation as time pro- 
gresses is a consequence of the subsequent scattering events, 
which leads to the correlation between paired modes becom- 
ing essentially zero at late times. 

Note that in Fig. 1 1 we have only plotted results for which 
Nj > 1/ 4. Although we can calculate the correlation function 
for all modes and at all times, the results have no physical 
meaning for unpopulated modes. Furthermore, modes with 
very low population produce very noisy results. 

The analogue to Eq. (76) for measuring the population cor- 
relation between paired modes is 



where for conciseness we have written fy-p = ^-p (x, t) and 
S-p = 8-p (x, x). As with ' (t), Eq. (75) should return unity 
for regions of coherently distributed density and two for re- 
gions of Gaussian distributed density. 

In Fig. 10 we plot (x, t) for the collision, calculated us- 
ing Eq. (75), at a sequence of times. These plots again show a 
uniformly coherent initial state, as required, and as time pro- 
gresses we observe that close to the centre of the collision 
the system becomes increasing incoherent, reflecting the cre- 
ation of the halo, so that by 37.7 ms the central region shows 
5 (2) (x,t)~2. 



(2) 

9ij 



(t) 



a\a\aiaj) (t) 



alai)(t)(a]a 3 )(t) 



(78) 



We have calculated this quantity, and have observed essen- 
tially the same behavior as for the amplitude correlations, 
i.e. a definite correlation exists at early times, that decreases 
rapidly with time. We note however that a very large error ex- 

(2) 

ists in the calculations of g\- (of order ±2 at t = 6 ms), much 
larger than for the amplitude correlations, so that quantifying 
the degree of this correlation is difficult. 
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FIG. 1 1 : Magnitude of the mode amplitude pair creation correlation 
function glj' 2 ^ (i) on the plane v z = at 6.6 ms (a), 13.1 ms (b), 
19.7 ms (c) and 37.7 ms (d). Results are only shown for those modes 
with Nj > 1 /4. 



E. Autocorrelation function 

As shown by Figs. 1, 2 and 3, the scattering halo in mo- 
mentum space is characterized by discrete phase grains. To 
quantify the size of these grains, and hence the range of phase 
order within the halo, we use the normalized autocorrelation 
function 



$ (t) 



a\a 1 ) (t) 
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which in terms of moments of the Wigner function is 
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Although modes i and j can be arbitrarily chosen, for this 
analysis we fix mode i at Vj = (0, 1.8 mms -1 , l) and vary j 
over the full mode space, for which the autocorrelation func- 
tion is shown on the planes v z = and v x = in Fig. 12. 
From these plots we observe that the significant feature of the 
autocorrelation function is a small ellipsoid, visible as a black 
grain centered on Vj = v;, whose size and shape remains 
relatively unchanged through time. 

To more precisely quantify the time-dependent shape of 
this grain we have fitted a three-dimensional Gaussian pro- 
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FIG. 12: Normalized mode space autocorrelation function g^ (t) 
on the planes v z = (a,c,e,g) and = (b,d,f,h) at 6.6 ms 
(a,b), 13.1 ms (c,d), 19.7 ms (e,f) and 37.7 ms (g,h). Here v; = 
(0, 1.8 mm s _1 , 0) and results are only shown for modes with Nj > 
1/4. 



file to the correlation function, over a region of approximately 
the same extent as the grain. At the time when significant 
population is first established in the halo (8 ms), the HWFM 
widths of the fitted Gaussian are approximately 0.3 mms -1 , 
0.3 mms -1 and 0.8 mms -1 in the v x , v y and v z directions 
respectively. These widths reflect the asymmetry of the initial 
condensate wavepackets (in velocity space), and are approxi- 
mately 1.5 times as large as the velocity radii of the conden- 
sate wavepackets in the Thomas-Fermi approximation (being 
0.19 mms -1 in the v x and v y directions and 0.54 mms -1 in 
the v z direction). At later times, the fitted widths are found 
to increase anisotropically, so that by 40 ms the fitted widths 
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FIG. 13: Total coherent and incoherent populations calculated us- 
ing the extrapolation method for varying collision parameters. Plots 
(a,b) show (respectively) the coherent and incoherent populations 
for collisions with Av = 4.0 mms -1 with (lowest to highest) 
iV (t = 0) = {1, 2, 3, 4} x 10 6 . Plots (c,d) show the coherent and 
incoherent populations for collisions with No (t = 0) = 2 x 10 6 
and (highest to lowest in (c) and lowest to highest in (d)) Aw = 
{1.0,2.0,3.0,4.0} rams" 1 . 



this peak occurs decreases, with increasing particle number. 
The dashed parts of the curves for N (t = 0) = {3, 4} x 10 6 
indicate results where the real particle density has begun to 
exit the simulation volume. 



H. Collision speed 

In Fig. 13 (c,d) we plot the total coherent and incoherent 
populations for collisions with No(t = 0) = 2 x 10 6 and 
Av = {1.0, 2.0, 3.0, 4.0} mm s _1 . These curves show that the 
rate of incoherent population scattering increases (again found 
to be roughly linearly) with increasing relative wavepacket 
speed, a result which is expected from classical collisional 
theory (nav). However, as can be seen for the curves with 
Av = {3.0, 4.0} mm s~\ the total amount of incoherent pop- 
ulation generated is less strongly dependent upon Av, as the 
overlap time of the condensate wavepackets, and hence the 
period over which the major mechanism of incoherent parti- 
cle formation acts, scales as « 1/yAv. For the lower values 
of Av, broadening of the condensate wavepackets leads to the 
population exiting the simulation volume before packet sepa- 
ration. Note that the the speed of sound for this system (taken 
as for a homogeneous system at the peak density of the con- 
densate at t = 0) is 1.5 mms -1 . 



are approximately 0.3 mms -1 , 0.6 mms -1 and 0.9 mms -1 . 
These increases are driven by scattering events, and reflect the 
inverse extent of the colliding system in coordinate space (see 
Fig. 4). 



F. Parameter dependence 

In this section we explore the dependence of the coherent 
and incoherent populations on the initial condensate number 
and on the initial relative wavepacket speed. These parame- 
ters, which are relatively easily adjusted experimentally, are 
the only ones we change — all other simulation parameters, 
including Uq, lu x . v ,z and the details of the simulation grids, 
are held fixed. We obtain the time-dependent total coherent 
and incoherent populations for each parameter set using the 
extrapolation method outlined in section VB3. For all pa- 
rameter sets two trajectories were simulated. 



G. Initial condensate number 

In Fig. 13 (a,b) we plot the total coherent and incoherent 
populations for collisions at Av = 4.0 mms -1 with total ini- 
tial condensate numbers of N (t = 0) = {1, 2, 3, 4} x 10 6 . 
From these curves we can see that by increasing the total 
number of particles within the system, the amount of inco- 
herent population generated by the collision also increases. 
We find that the peak rate of incoherent particle formation in- 
creases roughly linearly with No (0), and that the time that 



I. The halo formation process 

The basic mechanism of halo formation is pairwise scat- 
tering of condensate particles into two distinct halo modes, 
which can be viewed as a four-wave mixing process. How- 
ever, not all modes within the halo region exhibit population 
growth, a feature that can be understood in terms of the phase 
relationships required for growth. For plane wave modes, the 
nonlinear portion of the mode amplitude evolution is 



da (nonlin) jj 
I = > 



k r ,k s +k t , 



(81) 



rsteL 



representing a process where s + t — > j + r. Writing the mode 
amplitudes as otj (t) = \fJTj cxp [iOj (t)], the rate of popu- 
lation change for a single mode (in the absence of external 
potentials) is 



dhfn 2U 



dt 



V 



^/NjNrNsNt Sin (Qjrst) Skj +k r ,k s +k t : 



rsteL 



(82) 

where ®j rs t = Oj + 6 r — S — 8 t . The phase evolution corre- 
sponding to Eq. (82) is 



dt 




COS {Qjrst) Skj+K^s+kf 



(83) 

At early times, the evolution of the halo modes is driven al- 
most entirely by the condensate wavepackets, with the dom- 
inant processes being scattering events involving the destruc- 
tion of one quanta from each wavepacket. The role of the 
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phase Qj rs t is critical, as can be understood from a simplified 
model of four modes only, two highly occupied modes s and 
t representing the two condensate packets, and a pair of scat- 
tered modes j and r of lower occupation, for which k r = — k, 
(in the centre of mass frame). Initially the phase <dj rs t is ran- 
domly set, and if — tt < <dj rs t < 0, population transfers from 
the pair (s,t) into the pair (j,r), while if < Qjrst < tt, 
the population transfer is in the opposite direction. The phase 
Qj rs t evolves towards either — n/2 (maximum gain for scat- 
tered modes) or +tt/2 (maximum loss), where it stabilizes 
while population is still available for transfer. Thus some 
pairs (j, r) will grow, while others will decrease in size. If 
JVj <C M r , the phase 9j can change rapidly (see Eq. (83)), 
possibly leading to a change in the direction of population 
transfer. 

For the multimode case, the main additional feature is that 
the scattered modes no longer need be precise momentum op- 
posites, because of the range of momentum modes available in 
the condensate packets. Thus for a particular scattered mode 
j there is a range of possible modes near its conjugate mo- 
mentum mode r which can be convolved with the conden- 
sate packets to contribute to the growth (or loss) in j, (see 
Eq. (81)). Labelling this set of modes about r as R, then from 
the properties of convolution, the size of R is about twice that 
of a condensate packet. There is of course a corresponding set 
of modes J about j, and if the phases and populations for a 
pair of modes from R and J are favorable for growth, these 
phases can quickly be locked across the regions R and J (see 
Fig 3). Population growth is a stimulated process, so while 
R and J will now grow rapidly, regions which have not estab- 
lished a favorable phase are left behind. 

We note that this discussion helps explain the incomplete 
pair correlation observed for halo modes j and r, because the 
scattered pair are not required to have exactly opposite mo- 
mentum, due to the range of momentum modes available in 
the condensate packets. 

VI. CONCLUSION 

In this paper we have presented a complete derivation of 
the truncated Wigner method, paying particular attention to 
the limits of validity of the Wigner truncation. Using this 
formalism, we have presented simulation results, both from 
a single trajectory and from ensembles of trajectories, of col- 
lisions between distinct Bose-Einstein condensate wavepack- 



ets. This formalism includes both stimulated and spontaneous 
processes, allowing for a complete treatment of system dy- 
namics. In particular we have observed the generation of 
highly populated s-wave scattering haloes. Previous treat- 
ments of similar collisional systems have been restricted to 
the low-scattering limit, and have not provided a complete de- 
scription of the field. 

The most significant process not included in our treatment 
of these collisions is three-body recombination [37]. In a fur- 
ther paper [38] we extend the truncated Wigner method to in- 
clude this process, and investigate its effect on a simple sys- 
tem. For the colliding systems considered in this paper we 
have found using this extended formalism that three-body re- 
combination has a negligible effect on the dynamics. 



A. Relationship to the projected Gross-Pitaevskii equation 

The various forms of the truncated Wigner differential 
equation, Eqs. (36,37,38), are functionally identical to the 
Projected Gross-Pitaevskii equation of Davis et al. [16, 17, 
21]. However the projected Gross-Pitaevskii equation treat- 
ment has been used only for relatively high temperature situ- 
ations, in which the thermal fluctuations are much larger than 
the quantum fluctuations, whereas in the truncated Wigner 
method the quantum mechanical nature of the system is still 
largely present in the form of mode amplitude fluctuations in 
the initial state (see below). The initial state of the Gross- 
Pitaevskii equation system would represent the limit, in the 
truncated Wigner function method, in which the size of the 
quantum fluctuations becomes negligible in comparison with 
the size of the order parameter of the field, while keeping the 
number of particles in the field constant. Thus the Gross- 
Pitaevskii equation can be considered to represent the thermo- 
dynamic limit of the truncated Wigner approach where spon- 
taneous, and spontaneously initiated, processes are far less 
important than stimulated processes. Of course the alternate 
approach is then simply to take the Gross-Pitaevskii equa- 
tion, and "seed" each mode with a certain amount of noise, 
such that previously unavailable spontaneous processes be- 
come possible. While this approach would avoid some of the 
more complicated aspects of the truncated Wigner method, 
one could not then make the unambiguous connections to the 
underlying quantum nature of the system which the truncated 
Wigner method allows. 
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